%PricingInvestment
% Pricing regressions
% Liuren Wu, liurenwu@gmail.com
format compact; clearvars;
load('../data/Momentsw.mat');
load('../data/PLw.mat','PLD');

    tic;
%BMS greeks, with F normalized to 100
[Theta, dVega, dGamma, dVanna, dVolga] = FunForwardVGVVGreeks(100, d1w, Iw, matw);


d1m=repmat(d1(:)',nm,1);d1m=(d1m(:)); nd1m=normcdf(d1m);
Nd1m=normcdf(d1m);
mm=repmat(mat(:),1,nk);mm=(mm(:));smm=sqrt(mm);lmm=log(mm);

pdm=repmat(pdelta(:)'*100,nm,1); pdm=pdm(:);

nmk=nm*nk;%cross section size
bdsd=2;bdsm=1;

np=4;
R2=NaN(Tw,nh,nc);
BV=NaN(Tw,np,nh,nc);
B2V=NaN(Tw,np,nc);B2T=B2V;
RP=NaN(Tw,np,nc);RRP=RP;NRP=RP;
R22=NaN(Tw,nc);AR2=R22;
ev=NaN(Tw,nmk,nc);evv=ev;
CSMoments=NaN(Tw,nmk,4,nh+1,nc);
FPL=NaN(Tw,np+1,3,nc);
RPV=NaN(Tw,np,3,nc);%raw,risk, notional
VC=NaN(Tw,np+1,nc);
BO=NaN(Tw,nc);
NOBV=NaN(Tw,nh+1,nc);
FactorModel=NaN(Tw,np+2,nc);
RVC=NaN(Tw,np+1,nc);ARVC=RVC;
AFactorModel=FactorModel;
wr=NaN(nmk,np+1,Tw);
wx=NaN(nmk,np,Tw);
nz=4;
wz=NaN(nmk,nz,Tw);
Aev=ev;AB2V=NaN(Tw,nz,nc);
AVC=NaN(Tw,nz+1,nc);
RRR=NaN(Tw,2,nc);
RNDS=NaN(nmk,Tw);
IOX=zeros(nmk,Tw);
MSEOS=NaN(Tw,16,nc);
MSEOZ=NaN(Tw,16,nc);
AVGM=NaN(Tw,4,nc);
ARRP=NaN(Tw,np,nc);
tradeflag=NaN(Tw,nc);
AFPL=FPL;
rng('default');
for t=54:Tw
    %add random number
    rdt=rand(nmk,1);
    RNDS(:,t)=rdt;
    idxis=rdt<=prctile(rdt,80);
    idxos=idxis==0;
    IOX(:,t)=idxis;

    for c=1:nc
        y=reshape(-Theta(t,:,:,c),nmk,1);
        vegat=reshape(dVega(t,:,:,c),nmk,1);
        gammat=reshape(dGamma(t,:,:,c),nmk,1);
        vannat=reshape(dVanna(t,:,:,c),nmk,1);
        volgat=reshape(dVolga(t,:,:,c),nmk,1);
        It=reshape(Iw(t,:,:,c),nmk,1);

        xx=NaN(nmk,np,nh);
        mom=NaN(nmk,4,nh);
        rhov=NaN(nh,1);
        for h=1:nh
            hmt=reshape(HMw(t,:,:,c,h),nmk,1);
            hvt=reshape(HVw(t,:,:,c,h),nmk,1);
            hgt=reshape(HGw(t,:,:,c,h),nmk,1);
            hot=reshape(HOw(t,:,:,c,h),nmk,1);

             hmt= FunMomentCSSmoothing(hmt,Nd1m,lmm,bdsd,bdsm);
             hvt= FunMomentCSSmoothing(hvt,Nd1m,lmm,bdsd,bdsm);
             hot= FunMomentCSSmoothing(hot,Nd1m,lmm,bdsd,bdsm);
             hgt= FunMomentCSSmoothing(hgt,Nd1m,lmm,bdsd,bdsm);

            rhov(h)=nanmean(hgt(:))./sqrt(nanmean(hvt(:)).*nanmean(hot(:)));
            hovt=sqrt(hot.*hvt);

            mom(:,:,h)=[hmt,hvt,hgt,hot];
            X=[vegat.*(hot),  0.5*gammat.*hvt, 0.5*volgat.*hot,vannat.*(hovt)];
            ii=isfinite(sum(X,2)+y);
            Ni=sum(ii);
            NOBV(t,h,c)=Ni;
            if Ni>15
                yi=y(ii); Xi=X(ii,:);
                B=Xi\yi;
                ei=yi-Xi*B;  R2t=1-mean(ei.^2)/var(yi);

                R2(t,h,c)=R2t;
                BV(t,:,h,c)=B;
                xx(ii,:,h)=Xi;
            else
                %disp(datestr(wdate(t)))
                %Ni
                %pause
            end
        end
        [R2t,indmr]=max(R2(t,:,c)); w=zeros(3,1); w(indmr)=1;
        mom(:,:,4)=mom(:,:,1)*w(1)+mom(:,:,2)*w(2)+mom(:,:,3)*w(3);
        CSMoments(t,:,:,:,c)=mom;
        X=xx(:,:,1)*w(1)+xx(:,:,2)*w(2)+xx(:,:,3)*w(3);
        
        Bnull=[0,1,1,rhov(1)]';
        BO(t,c)=rhov(1);
        ii=isfinite(sum(X,2)+y);
        Ni=sum(ii);NOBV(t,h+1,c)=Ni;
        if Ni>15
           tradeflag(t,c)=1; 
           AVGM(t,:,c)=mean(mom(ii,:,4)); 
           
           iis=ii&idxis; ios=ii&idxos;
           yi=y(ii);
            Xi=X(ii,:);
            B=Xi\yi;
            e=y-X*B;ev=e./gammat./It;
            ei=e(ii); eiv=ev(ii); %full sample
            eios=e(ios);eivos=ev(ios); %full sample estimates on out of sample
           
            %%in-sample out of sample
            yii=y(iis);
            Xii=X(iis,:);
            Bi=Xii\yii;
            e=y-X*Bi;
            evs=e./gammat./It;
            eis=e(iis); eos=e(ios);
            evis=evs(iis); evos=evs(ios); 

            
            
            MSEOS(t,:,c)=[mean(ei.^2), mean(eis.^2), mean(eos.^2),mean(eios.^2), ...
                mean(yi.^2), mean(y(iis).^2), mean(y(ios).^2), mean(y(ios).^2), ...
                mean(eiv.^2),mean(evis.^2),mean(evos.^2),mean(eivos.^2), ...
                mean(abs(eiv)),mean(abs(evis)),mean(abs(evos)),mean(abs(eivos))];

            %%

            
            R2t=1-mean(ei.^2)/mean(yi.^2);
            AR2t=1-(sum(ei.^2)./(Ni-np))/var(yi);
            
            eiv=ei./gammat(ii)./It(ii);
            
            sb=sqrt(diag(inv(Xi'*Xi)).*var(ei));
            tb=B./sb;
            
            wx(ii,:,t,c)=Xi/(Xi'*Xi);
            cXX=(Xi'*Xi);cy=yi'*yi;
            aat= [((cXX*B/cy).*B);  R2t];
            VC(t,:,c)=aat;
            
            
            R22(t,c)=R2t;
            AR2(t,c)= AR2t;
            
            
            B2V(t,:,c)=B; B2T(t,:,c)=tb;
            ev(t,ii,c)=ei;
            evv(t,ii,c)=eiv;
            
            RRP(t,:,c)=B-Bnull;
            Hr=[Xi, ei];              portwr=Hr/(Hr'*Hr);
            sH=sqrt(mean(Hr.^2));     portwh=portwr.*sH;   RP(t,:,c)=RRP(t,:,c).*sH(1:np);
            swrt=sum(abs(portwr));    portwn=portwr./swrt; NRP(t,:,c)=RRP(t,:,c)./swrt(1:np);
            
            wr(ii,:,t,c)=portwr;
            
            RPV(t,:,:,c)=[RRP(t,:,c)',RP(t,:,c)',NRP(t,:,c)'];
            
            % pricing under common risk assumption
            Z=[vegat, 0.5*gammat, 0.5*volgat,vannat];
            yi=y(ii);
            Zi=Z(ii,:);
            C=Zi\yi;
            ez=y-Z*C; ezv=ez./gammat./It;
            eiz=ez(ii);eizv=ezv(ii);
            eizos=ez(ios);eizvos=ezv(ios); %full sample estimates on out of sample
           
            %out of sample
            yii=y(iis);
            Zii=Z(iis,:);
            Ci=Zii\yii;
            e=y-Z*Ci;
            evs=e./gammat./It;
            eis=e(iis); eos=e(ios);
            evis=evs(iis); evos=evs(ios);
          
            MSEOZ(t,:,c)=[mean(eiz.^2),mean(eis.^2), mean(eos.^2),mean(eizos.^2), ...
                mean(yi.^2), mean(y(iis).^2), mean(y(ios).^2), mean(y(ios).^2), ...
                mean(eizv.^2),mean(evis.^2),mean(evos.^2),mean(eizvos.^2), ...
                mean(abs(eizv)),mean(abs(evis)),mean(abs(evos)),mean(abs(eizvos))];


            
            R2j=1-mean(eiz.^2)/mean(yi.^2);
            cZZ=(Zi'*Zi);cy=yi'*yi;
            AVC(t,:,c)= [((cZZ*C/cy).*C);   R2j];
            
            ammt= mean(mom(ii,:,4)); 
            ammt=mom(13,:,4);
            CHNull=[0,ammt(2), ammt(4), ammt(3)]';
            ARRP(t,:,c)=C-CHNull; %risk-neutral v historical
            
            AB2V(t,:,c)=C;
            Aev(t,ii,c)=eiz;
            wz(ii,:,t,c)=Zi/(Zi'*Zi);
            
            %return factor model
            tt=t-52:t-1;
            vrt=reshape(nanstd(PLD(tt,:,:,c)),nmk,1);
            

            PLt=reshape(PLD(t,:,:,c),nmk,1);
            PLi=PLt(ii);
            ij=isfinite(PLi);
            if sum(ij)>15
                %regression
                Hi=[Xi,ei];
                Hj=Hi(ij,:);yj=-PLi(ij);
                Bj=Hj\yj;
                ej=yj-Hj*Bj;R2j=1-mean(ej.^2)/mean(yj.^2);
                
                FactorModel(t,:,c)=[Bj;R2j];
                
                cXX=(Hj'*Hj);cy=yj'*yj;
                aat= [((cXX*Bj/cy).*Bj)];
                RVC(t,:,c)=aat;
                
                
                for k=1:np+1 %risk-adjusted factor portfolio
                    w0=portwr(:,k);
                    w1=portwh(:,k);
                    w2=portwn(:,k);
                    ww=[w0,w1,w2];

                    FPL(t,k,:,c)=-PLi(ij)'*ww(ij,:); %PL for next day, recorded on today
                end
               
                % Hrr=[Xi, ei]./vrt(ii);              
                % portwwr=(Hrr/(Hrr'*Hrr))./vrt(ii);swwrt=sum(abs(portwwr));    
                % portwwn=portwwr./swwrt;
                % k=np+1;
                % FPL(t,k,3,c)=-PLi(ij)'*portwwn(ij,k);

                %return factor model under common risk assumption
                Hz=[Zi, eiz];
                Hzj=Hz(ij,:);yjj=-PLi(ij);
                Bz=Hzj\yjj;
                ezj=yjj-Hzj*Bz;R2z=1-mean(ezj.^2)/mean(yjj.^2);
                AFactorModel(t,:,c)=[Bz;R2z];
        
                cZZ=(Hz'*Hz);cy=yj'*yj;
                aat= [((cZZ*Bz/cy).*Bz)];
                ARVC(t,:,c)=aat;
                RRR(t,:,c)=[R2j,R2z];

                %factor portfolio z
                portwr=Hz/(Hz'*Hz);
                sH=sqrt(mean(Hz.^2));     portwh=portwr.*sH; 
                swrt=sum(abs(portwr));    portwn=portwr./swrt;
                for k=1:np+1 %z factor portfolio
                    w0=portwr(:,k);
                    w1=portwh(:,k);
                    w2=portwn(:,k);
                    ww=[w0,w1,w2];
                    AFPL(t,k,:,c)=-PLi(ij)'*ww(ij,:);
                end


                
            end
            
            
            
            
            
            
        end
        
    end
end
toc;

save('../data/PricingResults4f2.mat','R22','VC','FPL','RRP','RPV','RP','NRP','CSMoments','BO','ev','evv','FactorModel', 'RVC', ...
    'AVC','AB2V','Aev','wr','wx','wz', 'ARVC','AFactorModel', 'MSEOS', 'MSEOZ','AVGM','tradeflag', 'AFPL','ARRP', ...
    'Theta','dVega','dGamma','dVanna','dVolga', ...
    'np','nmk','-mat');
%ResultsSummary;


return


